{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from [this Jupyter notebook](https://nbviewer.jupyter.org/github/krischer/seismo_live/blob/master/notebooks/Computational%20Seismology/The Finite-Difference Method/fd_ac1d.ipynb) by Heiner Igel ([@heinerigel](https://github.com/heinerigel)), Lion Krischer ([@krischer](https://github.com/krischer)) and Taufiqurrahman ([@git-taufiqurrahman](https://github.com/git-taufiqurrahman)) which is a supplemenatry material to the book [Computational Seismology: A Practical Introduction](http://www.computational-seismology.org/),  additional modifications by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving the 1D acoustic wave equation by finite-differences\n",
    "\n",
    "In the previous lectures we derived acoustic and elastic approximations to the partial differential equations of motion, describing seismic wave propagation in 3D isotropic elastic media. For homogeneous acoustic media we derived analytical solutions in terms of Green's functions. Finally, we dealt with the discretization of material parameters and wavefields in continous media and the finite-difference approximation of partial derivatives. \n",
    "\n",
    "In this lesson we wrap up the results of all previous lectures by solving the 1D acoustic wave equation for a homogeneous medium using the finite-difference method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite difference solution\n",
    "\n",
    "As derived in [this](https://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/01_Analytical_solutions/3_Acoustic_medium.ipynb) and [this lecture](https://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/01_Analytical_solutions/4_2D_1D_elastic_acoustic_approx.ipynb), the acoustic wave equation in 1D with constant density is\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,t)}{\\partial t^2} \\ = \\ vp(x)^2 \\frac{\\partial^2 p(x,t)}{\\partial x^2} + s(x,t) \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "with pressure $p$, acoustic velocity $vp$, and source term $s$. Both second derivatives can be approximated by a 3-point difference formula. For example for the time derivative, we get:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,t)}{\\partial t^2} \\ \\approx \\ \\frac{p(x,t+dt) - 2 p(x,t) + p(x,t-dt)}{dt^2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "and equivalently for the spatial derivative: \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,t)}{\\partial x^2} \\ \\approx \\ \\frac{p(x+dx,t) - 2 p(x,t) + p(x-dx,t)}{dx^2}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Injecting these approximations into the wave equation allows us to formulate the pressure p(x) for the time step $t+dt$ (the future) as a function of the pressure at time $t$ (now) and $t-dt$ (the past). This is called an **explicit scheme** allowing the $extrapolation$  of the space-dependent field into the future only looking at the nearest neighbourhood.\n",
    "\n",
    "After discretization of the P-wave velocity and pressure wavefield at the discrete spatial grid points $i = 0, 1, 2, ..., nx$ and time steps $n = 0, 1, 2, ..., nt$, we can replace the time-dependent part (upper index time, lower index space) by\n",
    "\n",
    "\\begin{equation}\n",
    " \\frac{p_{i}^{n+1} - 2 p_{i}^n + p_{i}^{n-1}}{\\mathrm{d}t^2} \\ = \\ vp_{i}^2 \\biggl( \\frac{\\partial^2 p}{\\partial x^2}\\biggr) \\ + s_{i}^n \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Solving for $p_{i}^{n+1}$ leads to the extrapolation scheme:\n",
    "\n",
    "\\begin{equation}\n",
    "p_{i}^{n+1} \\ = \\ vp_i^2 \\mathrm{d}t^2 \\left( \\frac{\\partial^2 p}{\\partial x^2} \\right) + 2p_{i}^n - p_{i}^{n-1} + \\mathrm{d}t^2 s_{i}^n\n",
    "\\end{equation}\n",
    "\n",
    "The spatial derivatives are determined by \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 p(x,t)}{\\partial x^2} \\ \\approx \\ \\frac{p_{i+1}^{n} - 2 p_{i}^n + p_{i-1}^{n}}{\\mathrm{d}x^2} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Eq. (1) is the essential core of the FD modelling code. Because we derived analytical solutions for wave propagation in a homogeneous medium, we should test our first code implementation for a similar medium, by setting:\n",
    "\n",
    "\\begin{equation}\n",
    "vp_i = vp0 \\notag\n",
    "\\end{equation}\n",
    "\n",
    "at each spatial grid point $i = 0, 1, 2, ..., nx$, in order to compare the numerical with the analytical solution. For a complete description of the problem we also have to define initial and boundary conditions. The **initial condition** is \n",
    "\n",
    "\\begin{equation}\n",
    "p_{i}^0 = 0, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "so the modelling starts with zero pressure amplitude at each spatial grid point $i = 0, 1, 2, ..., nx$. As **boundary conditions**, we assume \n",
    "\n",
    "\\begin{align}\n",
    "p_{0}^n = 0, \\nonumber\\\\\n",
    "p_{nx}^n = 0, \\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "for all time steps n. This **Dirichlet boundary condition**, leads to artifical boundary reflections which would obviously not describe a homogeneous medium. For now, we simply extend the model, so that boundary reflections are not recorded at the receiver positions.\n",
    "\n",
    "Let's implement it ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries (PLEASE RUN THIS CODE FIRST!) \n",
    "# ----------------------------------------------\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "# Show Plot in The Notebook\n",
    "matplotlib.use(\"nbagg\")\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Sub-plot Configuration\n",
    "# ----------------------\n",
    "from matplotlib import gridspec \n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "\n",
    "         # number of grid points in x-direction\n",
    "         # grid point distance in x-direction\n",
    "         # wave speed in medium (m/s)\n",
    "         # source location in grid in x-direction\n",
    "         # receiver location in grid in x-direction\n",
    "         # maximum number of time steps\n",
    "         # time step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Source time function\n",
    "\n",
    "To excitate wave propagation in our 1D homogenous model, we use the first derivative of the Gaussian:\n",
    "\n",
    "\\begin{equation}\n",
    "s(t) = -2 (t-t_0) f_0^2 exp(-f_0^2 (t-t_0)^2) \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "as source time function at the discrete source position $isrc$. Where $t_0$ denotes a time shift and $f_0$ the dominant frequency of the source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Plot Source Time Function \n",
    "# -------------------------\n",
    "\n",
    " # dominant frequency of the source (Hz)\n",
    " # source time shift\n",
    "\n",
    "print('Source frequency =', f0, 'Hz')\n",
    "\n",
    "# Source time function (Gaussian)\n",
    "# -------------------------------\n",
    "\n",
    "\n",
    "# 1st derivative of a Gaussian\n",
    "src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "# Plot position configuration\n",
    "# ---------------------------\n",
    "plt.ion()\n",
    "fig1 = plt.figure(figsize=(6, 5))\n",
    "gs1  = gridspec.GridSpec(1, 2, width_ratios=[1, 1], hspace=0.3, wspace=0.3)\n",
    "\n",
    "# Plot source time function\n",
    "# -------------------------\n",
    "ax1  = plt.subplot(gs1[0])\n",
    "ax1.plot(time, src) # plot source time function\n",
    "ax1.set_title('Source Time Function')\n",
    "ax1.set_xlim(time[0], time[-1])\n",
    "ax1.set_xlabel('Time (s)')\n",
    "ax1.set_ylabel('Amplitude')\n",
    "\n",
    "# Plot source spectrum\n",
    "# --------------------\n",
    "ax2  = plt.subplot(gs1[1])\n",
    "      # source time function in frequency domain\n",
    "      # time domain to frequency domain\n",
    "      # plot frequency and amplitude\n",
    "      # only display frequency from 0 to 250 Hz\n",
    "ax2.set_title('Source Spectrum')\n",
    "ax2.set_xlabel('Frequency (Hz)')\n",
    "ax2.set_ylabel('Amplitude')\n",
    "\n",
    "ax2.yaxis.tick_right()\n",
    "ax2.yaxis.set_label_position(\"right\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytical Solution\n",
    "\n",
    "In [this lecture](https://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/01_Analytical_solutions/5_Greens_function_acoustic_1-3D.ipynb) we calculated the Green's functions for the homogenous acoustic wave equation\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2}{\\partial t^2} G(x,t;x_s, t_s) \\ - \\ vp_0^2 \\Delta G(x,t;x_s, t_s) \\ = \\delta (x-x_s) \\delta (t-t_s) \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "where $\\Delta$ denotes the Laplace operator and the $\\delta-$function is defined as\n",
    "\n",
    "\\begin{equation}\n",
    "\\delta(x) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "\t\\infty &x=0 \\\\\n",
    "\t0 &x\\neq 0 \t\n",
    "\\end{array}\n",
    "\\right.\\notag\n",
    "\\end{equation}\n",
    "\n",
    "with the constraint:\n",
    "\n",
    "\\begin{equation}\n",
    "\\int_{-\\infty}^{\\infty} \\delta(x)\\; dx = 1.\\notag\n",
    "\\end{equation}\n",
    "\n",
    "When comparing numerical with analytical solutions the functions that - in the limit - lead to the $\\delta-$function will become very important. An example is the boxcar function\n",
    "\n",
    "$$\n",
    "\\delta_{bc}(x) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "\t1/dx &|x|\\leq dx/2 \\\\\n",
    "\t0 &\\text{elsewhere} \t\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "fulfilling these properties as $dx\\rightarrow0$. These functions are used to properly scale the source terms to obtain correct absolute amplitudes.\n",
    "\n",
    "To describe analytical solutions for the acoustic wave equation we also make use of the unit step function, also known as the Heaviside function, defined as \n",
    "\n",
    "$$\n",
    "H(x) = \\left\\{\n",
    "\\begin{array}{ll}\n",
    "    0 &x<0 \\\\\n",
    "\t1 &x \\geq 0 \t\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "The Heaviside function is the integral of the $\\delta-$function (and vice-versa the $\\delta$-function is defined as the derivative of the Heaviside function). In 1D case, the Greens function is proportional to a Heaviside function. \n",
    "\n",
    "$$\n",
    "G=\\frac{1}{2vp_0}H\\left(t-\\frac{|x|}{vp_0}\\right)\n",
    "$$\n",
    "\n",
    "As the response to an arbitrary source time function can be obtained by convolution \n",
    "\n",
    "$$\n",
    "G_{seis} = G(x,t;x_s,t_s) * s(t),\n",
    "$$\n",
    "\n",
    "this implies that the propagating waveform is the integral of the source time function. The response is shown for a source time function with a 1st derivative of a Gaussian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Analytical solution\n",
    "# -------------------\n",
    "G    = time * 0.\n",
    "\n",
    "# Initialize coordinates\n",
    "# ----------------------\n",
    "x    = np.arange(nx)\n",
    "x    = x * dx       # coordinate in x-direction\n",
    "\n",
    "for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "    if (time[it] - np.abs(x[ir] - x[isrc]) / vp0) >= 0:\n",
    "        G[it] = 1. / (2 * vp0)\n",
    "        \n",
    "    # get limit value from the maximum amplitude\n",
    "\n",
    "# Plotting convolution of Green's function with source wavelet\n",
    "plt.plot(time, Gc)\n",
    "plt.title(\"Analytical solution\" )\n",
    "plt.xlabel(\"Time, s\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of numerical with analytical solution\n",
    "\n",
    "In the code below we solve the homogeneous acoustic wave equation by the 3-point difference operator and compare the numerical results with the analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Plot Snapshot & Seismogram (PLEASE RERUN THIS CODE AGAIN AFTER SIMULATION!) \n",
    "# ---------------------------------------------------------------------------\n",
    "\n",
    "# Initialize empty pressure arrays\n",
    "# --------------------------------\n",
    " # p at time n (now)\n",
    " # p at time n-1 (past)\n",
    " # p at time n+1 (present)\n",
    " # 2nd spatial derivative of p\n",
    "\n",
    "# Initialize model (assume homogeneous model)\n",
    "# -------------------------------------------\n",
    "\n",
    "      # initialize wave velocity in model\n",
    "\n",
    "# Initialize empty seismogram\n",
    "# ---------------------------\n",
    " \n",
    "\n",
    "# Plot position configuration\n",
    "# ---------------------------\n",
    "plt.ion()\n",
    "fig2  = plt.figure(figsize=(6, 4))\n",
    "gs2   = gridspec.GridSpec(1, 2, width_ratios=[1, 1], hspace=0.3, wspace=0.3)\n",
    "\n",
    "# Plot 1D wave propagation\n",
    "# ------------------------\n",
    "# Note: comma is needed to update the variable\n",
    "ax3  = plt.subplot(gs2[0])\n",
    "leg1,= ax3.plot(isrc, 0, 'r*', markersize=11) # plot position of the source in snapshot\n",
    "leg2,= ax3.plot(ir, 0, 'k^', markersize=8) # plot position of the receiver in snapshot\n",
    "up31,= ax3.plot(p) # plot pressure update each time step\n",
    "ax3.set_xlim(0, nx)\n",
    "ax3.set_ylim(-lim, lim)\n",
    "ax3.set_title('Time Step (nt) = 0')\n",
    "ax3.set_xlabel('nx')\n",
    "ax3.set_ylabel('Amplitude')\n",
    "ax3.legend((leg1, leg2), ('Source', 'Receiver'), loc='upper right', fontsize=10, numpoints=1)\n",
    "\n",
    "# Plot seismogram \n",
    "# ---------------\n",
    "# Note: comma is needed to update the variable\n",
    "ax4  = plt.subplot(gs2[1])\n",
    "leg3,= ax4.plot(0,0,'r--',markersize=1) # plot analytical solution marker\n",
    "leg4,= ax4.plot(0,0,'b-',markersize=1) # plot numerical solution marker\n",
    "up41,= ax4.plot(time, seis) # update recorded seismogram each time step\n",
    "up42,= ax4.plot([0], [0], 'r|', markersize=15) # update time step position\n",
    "ax4.yaxis.tick_right()\n",
    "ax4.yaxis.set_label_position(\"right\")\n",
    "ax4.set_xlim(time[0], time[-1])\n",
    "ax4.set_title('Seismogram')\n",
    "ax4.set_xlabel('Time (s)')\n",
    "ax4.set_ylabel('Amplitude')\n",
    "ax4.legend((leg3, leg4), ('Analytical', 'FD'), loc='upper right', fontsize=10, numpoints=1)\n",
    "\n",
    "plt.plot(time,Gc,'r--') # plot analytical solution\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     42
    ]
   },
   "outputs": [],
   "source": [
    "# 1D Wave Propagation (Finite Difference Solution) \n",
    "# ------------------------------------------------\n",
    "\n",
    "# Calculate Partial Derivatives\n",
    "# -----------------------------\n",
    "for it in range(nt):\n",
    "    \n",
    "    # FD approximation of spatial derivative by 3 point operator\n",
    "    # ----------------------------------------------------------\n",
    "\n",
    "\n",
    "    # Explicit Time Extrapolation by 3 point operator\n",
    "    # -----------------------------------------------\n",
    "\n",
    "    # Add Source Term at isrc\n",
    "    # -----------------------\n",
    "    # Absolute pressure w.r.t analytical solution\n",
    "\n",
    "                \n",
    "    # Remap Time Levels\n",
    "    # -----------------\n",
    "    \n",
    "    # Output Seismogram\n",
    "    # -----------------\n",
    "    \n",
    "    # Update Data for Wave Propagation Plot\n",
    "    # -------------------------------------\n",
    "    idisp = 2 # display frequency\n",
    "    if (it % idisp) == 0:\n",
    "        ax3.set_title('Time Step (nt) = %d' % it)\n",
    "        up31.set_ydata(p)\n",
    "        up41.set_ydata(seis)\n",
    "        up42.set_data(time[it], seis[it])\n",
    "        plt.gcf().canvas.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* implementation of the 1D acoustic wave equation \n",
    "* understanding the input parameters for the simulation and the plots that are generated\n",
    "* modifying source and receiver locations and observing the effects on the seismograms\n",
    "* allowing you to explore the finite-difference method"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
